function [ys,params,check] = Tables2and3_steadystate(ys,exo,M_,options_)
% Inputs: 
%   - ys        [vector] vector of initial values for the steady state of
%                   the endogenous variables
%   - exo       [vector] vector of values for the exogenous variables
%   - M_        [structure] Dynare model structure
%   - options   [structure] Dynare options structure
%
% Output: 
%   - ys        [vector] vector of steady state values for the the endogenous variables
%   - params    [vector] vector of parameter values
%   - check     [scalar] set to 0 if steady state computation worked and to
%                    1 of not (allows to impose restrictions on parameters)
%
% Copyright (C) 2014-2020 Johannes Pfeifer
% 
%  This is free software: you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation, either version 3 of the License, or
%  (at your option) any later version.
% 
%  It is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
% 
%  For a copy of the GNU General Public License,
%  see <http://www.gnu.org/licenses/>.

% read out parameters to access them with their name

NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = M_.param_names{ii};
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
% initialize indicator
check = 0;


%% Enter model equations here
% The steady state values are derived from solving the steady state
% equilibrium of the model for both concentrated and low-concentration
% sectors, including the off-equilibrium path Bayesian-Nash equilibrium for
% the concentrated sector. We use the standard unit normalization for
% steady state productivitry levels and, which affects the scaling of the enbdogenous variables. 

%Steady State Values

 A1=1;
 A2=1;
K1=7.45e-5;
I1= delta1*K1;
K2=1.0949e-5;
I2=delta2*K2;
p1=0.099;
p2=0.062;
u1=0.807;
D1=(p1-h1)*(A1*(u1*K1)^psik1)-delta1*K1-m1*K1;
D2=(p2-h2)*(A2*K2^psik2)-delta2*K2-m2*K2;
S1=99*D1;
S2=99*D2;
EOmega=alph;
W=N1*D1+N2*D2;
P=(phi1^sigmaes*p1^(1-sigmaes)+phi2^sigmaes*p2^(1-sigmaes))^(1/(1-sigmaes));
c1=(W/P)*(P*phi1/p1)^sigmaes;
c2=(W/P)*(P*phi2/p2)^sigmaes;
C=(phi1*c1^((sigmaes-1)/sigmaes)+phi2*c2^((sigmaes-1)/sigmaes))^(sigmaes/(sigmaes-1));
U=C;
Q=U^(1-gamm);
I1tilde=delta1*7.29e-6;
S1tilde=99*3.29e-4;
R1=1/alph;
R2=1/alph;
Rf=1/alph;
Rf_real=Rf;
rf=Rf-1;
rf_real=rf;
Er1=R1-1;
Er2=R2-1;
Erp1=Er1-rf;
Erp2=Er2-rf;
c=log(C);
gc=0;
pi=log(P);
gpi=0;
IR1=I1/K1;
ir1=log(IR1);
gir1=0;
IR2=I2/K2;
ir2=log(IR2);
gir2=0;
gy1=0;
gy2=0;
w=log(W);
gw=0;
Y1=p1*A1*(u1*K1)^psik1;
Y2=p2*A2*K2^psik2;


       

%% end own model equations

params=NaN(NumberOfParameters,1);
for iter = 1:length(M_.params) %update parameters set in the file
  eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

NumberOfEndogenousVariables = M_.orig_endo_nbr; %auxiliary variables are set automatically
for ii = 1:NumberOfEndogenousVariables
  varname = M_.endo_names{ii};
  eval(['ys(' int2str(ii) ') = ' varname ';']);
end



